program fake_mad_gui
USE run_madx
USE pointer_lattice
implicit none
character*48 command_gino
type(layout), pointer :: R
TYPE(FIBRE),POINTER :: P 
INTEGER I,POS,k,j,mf,ks,no
LOGICAL(LP) DONEIT,orb,COSLIKE
REAL(DP) xv,js,jb,xs(2),xb(2),xt(lnv),xtn(lnv),d(3),angle(3)
type(pol_block) qf(2),qd(2),sexb
type(real_8) y(6)
real(dp) fix1(6),fix2(6),betax,closed(6),x(6),ent(3,3)
type(normalform) n1 
type(damap) id,nc,n,m,nr,A_t
type(pbfield) hrot
type(internal_state), target :: state
type(gmap) g,g0
integer PATCH_NEEDED
type(onelieexponent) h
type(pbresonance) hr
real(dp) prec
!real(dp) small_ray(6),small_number,mat(6,6),mat1(6,6),mat_exact(6,6),basis(3,3),omega(3)
!real(dp) lambda,x0(6),x1(6),x2(6),norm,x3(6),x4(6),mat3(6,6)
type(work) werk
character(nlp) sexname(100)
integer possex(100)
real(dp) ksex(100),ksex0(100)
type(c_damap) cmap,l
type(c_vector_field) f
type(c_normal_form) cnorm
type(damapspin) maps,ids
type(c_taylor) t
type(probe_8) ray8
type(probe) ray
character(5) c5
real(dp), allocatable :: ma(:,:),ve(:),gt(:,:)
integer neq,iter,kn,runs,niter


lmax=100.d0
use_info=.true.
my_estate=>state
prec=1.d-8 


!! This program reads a flat file provided by Chris Mayes
!! and produce by PTC via BMAD

! It runs the iterations in four groups 

! Runs=1 Fits in one iteration the tunes and all sextupole geometric resonances
! Runs=2 Fits first most deadly quartic resonances : interactive
! Runs=3 Fits additional quartic resonances
! Runs=4 Fits all the rest as well (difference resonances and tune shifts)

! The flat resulting file is printed "chrisflat_done.txt" and can be read by PTC; especially graphical interface is useful

! The file "ksex.txt" contains the Lie polynomial in phasors for easier checking that all is indeed zero
! It also contains the old and new sextupole strengths in PTC units

!!!! user stuff : one layout necessary before starting gui


 call ptc_ini_no_append


call READ_AND_APPEND_VIRGIN_general(M_U,'flat.txt')


c_verbose=.false.
runs=1
call kanalnummer(mf,"eq.txt")

r=>m_u%start

state=delta0

use_complex_in_ptc=my_true
call init(state,4,1)

call alloc(ids);call alloc(maps);
call alloc(ray8); call alloc(cmap);
call alloc(cnorm); call alloc(l); call alloc(f);
call alloc(t);

k=0
p=>r%start
do i=1,r%n
if(P%mag%p%nmul>2) then
if(p%mag%bn(3)/=0.d0) then
  write(6,*) i,p%mag%name,p%mag%bn(3)


 k=k+1
sexname(k)=p%mag%name
possex(k)=i 
ksex0(k)=p%mag%bn(3)
call add(p,3,0,0.d0)

endif
endif
p=>p%next
enddo

write(6,*) k, " sextupoles "

kn=k

2 continue 

neq= 12 
!niter=1  
if(runs==2) then 
 neq=18
endif
if(runs==3) then 
neq=24
endif
if(runs==4) then 
neq=31
endif
allocate(ma(neq,kn))
allocate(ve(neq+kn))
allocate(gt(neq+kn,neq+kn))

!do iter=1,niter

1 continue

ma=0.0_dp
ve=0.0_dp
gt=0.0_dp

do ks=1,kn
sexb=0
sexb%ibn(3)=1
sexb%name=sexname(ks)
r=sexb
 x=0.d0
 ray=0
 ray=x
 ids=1
 ray8=ray+ids

call track_probe(r,ray8,+state,fibre1=1)

 
maps=ray8
cmap=maps
call c_normal(cmap,cnorm,dospin=my_false,no_used=1)
!call print(cnorm%a_t,6)

cmap=cnorm%a_t**(-1)*cmap*cnorm%a_t
cmap=from_phasor(-1)*cmap*from_phasor(1)
call c_factor_map(cmap,l,f,1)
call c_clean_vector_field(f,f,prec)

 
c5='10001'
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)
ma(1,ks)=aimag(t).sub.'1'
ve(1)=aimag(t).sub.'0'

c5='00101'
t=f%v(3).par.c5
t=t<=5
if(c_verbose) call print(t,mf)
ma(2,ks)=aimag(t).sub.'1'
ve(2)=aimag(t).sub.'0'

c5='20000'     ! 1 0  J_x
t=f%v(1).par.c5 
t=t<=5
if(c_verbose) call print(t,mf)
ma(3,ks)=real(t).sub.'1'
ma(4,ks)=aimag(t).sub.'1'
ve(3)=real(t).sub.'0'
ve(4)=aimag(t).sub.'0'

c5='02000'     !  3 0
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)
ma(5,ks)=real(t).sub.'1'
ma(6,ks)=aimag(t).sub.'1'
ve(5)=real(t).sub.'0'
ve(6)=aimag(t).sub.'0'

c5='00200'     ! -1 2
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)
ma(7,ks)=real(t).sub.'1'
ma(8,ks)=aimag(t).sub.'1'
ve(7)=real(t).sub.'0'
ve(8)=aimag(t).sub.'0'

c5='00110'    !  1 0 J_y
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)
ma(9,ks)=real(t).sub.'1'
ma(10,ks)=aimag(t).sub.'1'
ve(9)=real(t).sub.'0'
ve(10)=aimag(t).sub.'0'

c5='00020'    ! 1 2
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(11,ks)=real(t).sub.'1'
ma(12,ks)=aimag(t).sub.'1'
ve(11)=real(t).sub.'0'
ve(12)=aimag(t).sub.'0'

if(neq==12) goto 100

c5='03000'    ! 4 0
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(13,ks)=real(t).sub.'1'
ma(14,ks)=aimag(t).sub.'1'
ve(13)=real(t).sub.'0'
ve(14)=aimag(t).sub.'0'

c5='01020'    ! 2 2
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(15,ks)=real(t).sub.'1'
ma(16,ks)=aimag(t).sub.'1'
ve(15)=real(t).sub.'0'
ve(16)=aimag(t).sub.'0'

c5='00030'    ! 0 4
t=f%v(3).par.c5
t=t<=5
if(c_verbose) call print(t,mf)


ma(17,ks)=real(t).sub.'1'
ma(18,ks)=aimag(t).sub.'1'
ve(17)=real(t).sub.'0'
ve(18)=aimag(t).sub.'0'

if(neq==18) goto 100

c5='12000'    ! 2 0 J_x
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(19,ks)=real(t).sub.'1'
ma(20,ks)=aimag(t).sub.'1'
ve(19)=real(t).sub.'0'
ve(20)=aimag(t).sub.'0'


c5='10020'    ! 0 2 J_x
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(21,ks)=real(t).sub.'1'
ma(22,ks)=aimag(t).sub.'1'
ve(21)=real(t).sub.'0'
ve(22)=aimag(t).sub.'0'


c5='01110'    ! 2 0  J_y
t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(23,ks)=real(t).sub.'1'
ma(24,ks)=aimag(t).sub.'1'
ve(23)=real(t).sub.'0'
ve(24)=aimag(t).sub.'0'

if(neq==24) goto 100

c5='21000'    ! J_x**2 tune 

t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(25,ks)=aimag(t).sub.'1'
ve(25)=aimag(t).sub.'0'


c5='10110'    ! J_x*J_y tune 

t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(26,ks)=aimag(t).sub.'1'
ve(26)=aimag(t).sub.'0'

c5='00210'    ! J_y**2 tune 

t=f%v(3).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(27,ks)=aimag(t).sub.'1'
ve(27)=aimag(t).sub.'0'


c5='01200'    ! 2 -2 

t=f%v(1).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(28,ks)=real(t).sub.'1'
ma(29,ks)=aimag(t).sub.'1'
ve(28)=real(t).sub.'0'
ve(29)=aimag(t).sub.'0'

c5='00120'    ! 0 2 J_y
t=f%v(3).par.c5
t=t<=5
if(c_verbose) call print(t,mf)

ma(30,ks)=real(t).sub.'1'
ma(31,ks)=aimag(t).sub.'1'
ve(30)=real(t).sub.'0'
ve(31)=aimag(t).sub.'0'

100 continue


write(6,*) cnorm%tune(1:2)
call kill_para(r)



write(6,*) "ks ",ks

!if(ks==1) goto 999
enddo


do i=1,kn
 gt(neq+i,neq+i)=1.0_dp
enddo

do i=1,neq
do ks=1,kn
 gt(i,neq+ks)=ma(i,ks)
 gt(neq+ks,i)=ma(i,ks)
enddo
enddo

close(mf)

xv=0
do i=1,12
xv=xv+abs(ve(i))
enddo
write(6,*) " Norm sex ",xv
xv=0
do i=13,neq
xv=xv+abs(ve(i))
enddo
write(6,*) " Norm oct ",xv




ks=neq+kn
call matinv(gt,gt,ks,ks,i)
write(6,*) " error = ",i

ve=-ve

ve=matmul(gt,ve)



k=1
p=>r%start
do i=1,r%n
 
IF(possex(k)==i) then
  call add(p,3,1,VE(neq+k))
  ksex(k)=p%mag%bn(3)
  k=k+1
endif
 
p=>p%next
enddo

if(runs>1) then
 write(6,*) " neq = ",neq
 write(6,*) " for more enter 1"
 read(5,*) iter
 if(iter==1) goto 1
endif
! enddo ! iter

deallocate(ma)
deallocate(ve)
deallocate(gt)

runs=runs+1
if(runs>4) goto 998
goto 2
998 continue
!goto 999
CALL  print_COMPLEX_SINGLE_STRUCTURE(r,"chrisflat_done.txt")

call kanalnummer(mf,"ksex.txt")

t=getpb(f,-1)
t=t/(-2*i_)
write(mf,*) " %%%%%%%%%%%%%%%%%%%%    Lie polynomial     %%%%%%%%%%%%%%%%%%%%% "
call print(t,mf)
write(mf,*) " %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% "



do i=1,kn
 write(MF,*) I,KSEX0(I),KSEX(I)
enddo

close(mf)

999 command_gino="OPENGINO"    
  call context(command_gino)   ! context makes them capital
  call call_gino(command_gino)  
 111 command_gino="MINI" 
  call context(command_gino)   ! context makes them capital
  call call_gino(command_gino)  
 



!!!!!!!!! vaguelynecessary baloney
command_gino="CLOSEGINO"
  call call_gino(command_gino)  



  call ptc_end
end program fake_mad_gui



